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NUMERICAL SIMULATION OF SOLAR CORONAL MAGNETIC FIELDS 


I. INTRODUCTION 

It is widely believed that the solar magnetic field is the driver for essentially all man- 
ifestations of coronal activity. Both the structure and the dynamics of coronal plasma are 
determined by the magnetic field. The basic idea is that photospheric footpoint motions 
stress the coronal field lines, thereby determining the coronal topology and producing mag- 
netic free energy that can be converted into plasma kinetic and thermal energy. Hence, the 
response of the coronal magnetoplasma to footpoint motions is of fundamental importance 
for understanding solar activity and many other astrophysical phenomena. 

In recent years there has been a great deal of work on this topic of coronal response to 
photospheric motions. Due to the complexity of the relevant equations, most of the effort 
has been on numerical simulations. We present in this paper the results of a numerical 
simulation that applies to the situation of twisting a flux tube in a magnetic arcade. The 
calculation is a dynamic, fully-three-dimensional MHD simulation in a geometry appropri- 
ate to the corona. To our knowledge this is the first one of its kind. The model and the 
numerical methods are described in the next Section. 


Manuscript approved August 2, 1990. 
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The results of our simulation are pertinent to several outstanding questions in solar 
coronal studies. One question concerns the topology of prominence magnetic fields. It is 
commonly believed that the cool, dense, H a emitting material in high quiescent promi- 
nences must be supported by the magnetic field (e.g. Tandberg-Hanssen 1974). In that 
case the magnetic field lines must form a “hammock” shape in the corona so that the 
cool material can rest at the bottom of the hammock. Several hammock-type models have 
been proposed (e.g. Kippenhahn and Schliiter 1957, Kuperus and Raadu 1974). The key 
question is, how does the field obtain such a geometry? A single current-free arcade has 
field lines that are concave-down in the corona. Perhaps the simplest and most appealing 
explanation is that photospheric motions twist the field lines of an initially current-free 
flux tube into the appropriate shape. Parker (1979) has shown that the twist will tend to 
propagate to the weakest part of the field, near the apexes of the flux tube. If sufficient 
twist is introduced, (of order one complete rotation or more), and if the twist is sufficiently 
localised at the apex, then some of the field lines must acquire a concave-up geometry. 
Analytic models for such a twisted flux tube have been constructed by several authors 
(e.g. Demoulin and Priest 1989). However, the analytic models consider only an equilib- 
rium flux tube in isolation and only for linear departures from a 2-dimensional equilibrium. 
They do not consider the effect of a realistic arcade structure for the field, the dynamic 
formation of the twists and the effect of instabilities. Our simulation allows us to study 
all these topological effects. 
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The role of ideal instabilities is an important issue in models for coronal activity. 
Since many observed phenomena such as prominence eruption, sprays and coronal mass 
ejections exhibit mass motions with velocities of the order of the coronal Alfven speed, it 
seems likely that they involve an ideal instability such as a kink. The basic idea is that if 
the coronal field is stressed sufficiently by photospheric motions, it becomes unstable and 
erupts outward carrying coronal and chromospheric plasma with it (e.g. Sturrock 1989). 
Two recent 3-d simulations seem to support this scenario (Strauss and Otani 1988; Mikic 
1990). These authors found that continuous twisting of a flux tube eventually resulted 
in a kink instability. However, their simulations were performed for a geometry in which 
the field is confined between two vertical plates and the system has cylindrical symmetry 
(before kinking). These constraints permit the buildup of arbitrarily large stresses . In the 
true coronal geometry the field is free to expand upward and there is no symmetry to the 
field. Hence the coronal magnetic field can relax gradually in response to photospheric 
motions. It is not clear if kinking will occur in this case as well, hence, more simulations 
are needed to determine whether ideal instabilities occur in more realistic coronal magnetic 
fields. 

Another issue is the role of magnetic reconnection. A number of authors (e.g. Pneu- 
man 1983, Moore 1988, Sturrock 1989, Van Ballegooijen and Martens 1989, etc) have 
invoked reconnection in a sheared arcade as a mechanism for both the formation and 
eruption of prominences. The concept underlying their models is that reconnection can 
redistribute the shear in an arcade such that part of the arcade flux becomes less sheared 
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while another part becomes progressively more sheared. In their model reconnection has 
the effect of concentrating the shear in a certain flux region. It is this highly sheared flux 
that eventually erupts, and leads to flares (Moore 1988). This is a very promising model. 
Simulations on a reduced two-dimensional model, coupled with the assumption that the 
reconnection did rearrange the shear as postulated, do show prominence eruption (Van 
Ballegooijen and Martens 1989). The full three-dimensional model has yet to be verified 
by rigorous calculations. However, there are serious doubts as to whether true reconnec- 
tion (as opposed to simple diffusion) can ever take place in a coronal arcade. Linear theory 
(e.g., Velli and Hood 1986) indicates that resistive tearing is suppressed by line-tying in a 
field that does not have a reversal of the unsheared initial field, as is the case for a dipole 
arcade. In the general non-linear regime, Greene (1988) has presented arguments for the 
absence of reconnection in a field with the topology of a single arcade. Hence, we believe 
that the role of reconnection in a coronal arcade is an open question which detailed, 3-d 
simulations can help to resolve. 

A related issue concerns the formation of DC electric qurrent sheets in closed magnetic 
configurations in the solar corona. Current sheets are very important for models of coronal 
heating. In order to account for the observed energy losses, it is widely believed that the 
electric currents in the corona must attain a highly filamented structure as in a current 
sheet. Parker (1972, 1987) proposes that current sheets form as a natural consequence 
of footpoint motions. He argues that a general 3-dimensional magnetic field such as that 
in the corona can have no well-behaved equilibrium state, and that magnetic singularities 
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corresponding to current sheets must form to allow the magnetic field to relax to a lower 
energy state. Several authors have numerically tested this hypothesis using 3-dimensional 
force-free-field codes (Van Ballegooijen 1988, Mikic et al, 1989). They concluded that 
their calculations did not support Parker’s hypothesis. These calculations, however, were 
performed for an initially uniform field confined between two plates. Our arcade geometry 
is more appropriate for the corona, and we include the dynamics, so that we can examine 
different aspects of the current-sheet formation question than the equilibria codes. 

In the next section we describe the exact model that we use for the simulation and 
our numerical code. The final section discusses the results and our conclusions. 

II. THE SIMULATION MODEL 

a) The Numerical Code 

We begin with the dissipative, incompressible magnetohydrodynamic (MHD) equa- 
tions, written in a dimensionless rotation form: 

| = vx«-vn+jxB + jjjV 2 v, ( la ) 

V • r = 0, (16) 

^ = vxB + V$-4vxVxA, (lc) 

at & 

and, 

V • B = 0, (Id) 
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where v(x,<) = (u,v,w) s flow velocity, w(x,<) = Vxv = vorticity, A(x,t) = magnetic 
vector potential, B (x,t) = VxA = magnetic field, j(x,i) = VxB = electric current 
density, II(x,i) = mechanical pressure + kinetic energy density, $(x,t) = magnetic scalar 
potential, S = Lundquist number, and M = viscous Lundquist number. Lundquist num- 
bers are Reynolds-like numbers, with the Alfven speed used for the characteristic velocity. 
We choose a gauge with $ = 0. In this representation, B is measured in terms of a char- 
acteristic field strength, e.g., B 0 = (v -1 (| B | 2 )) 2 , where the brackets <> indicate an 
integral over the system and V denotes the volume of integration. The velocities are mea- 
sured in units of the Alfven speed, Ca = Bo / (Air p) 3 where the mass density, p, is assumed 
to be constant and uniform. For a characteristic distance, Lo, the time is measured in 
units of the Alfven transit time, Lq/Ca • The Lundquist number S = CaLo/tj, where r? 
is the magnetic diffusivity. The viscous Lundquist number M — CaLo/v, where v is the 
kinematic viscosity. 

For the simulation reported here, 77 = .0025 and v — .2 so that S = 500 and M — 5. 
The low viscous Lundquist number was chosen to damp out fluid motions. In addition, 
it helps to dissipate small-scale magnetic structure via the Alfven effect, allowing us to 
perform fully 3D numerical simulations at a relatively large magnetic Lundquist number. 
However, such small values for the Lundquist numbers might be thought to preclude mag- 
netic reconnection a priori. To establish the possibility of reconnection at these Lundquist 
numbers, we have examined several linear and nonlinear problems where we have previ- 
ously observed magnetic reconnection. One such case is the linear magnetic reconnection 
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problem, which reduces to the MHD version of the Orr-Sommerfeld equation (Dahlburg 
et al., 1983, 1986). For the mean magnetic field B x (z) = tanh(z), we find that linearly 
unstable modes do exist at these Lundquist numbers. In terms of the wavenumber parallel 
to the mean magnetic field, a, and the growth rate, a;,-, we find, for example; (a,w,) = 
(0.6, 0.01), (1.0, 0.028057), (2.0, 0.03878), (3.0, 0.02255), and (4.0, 0.0096). This implies 
that linear instabilities are at least possible at these Lundquist numbers. We have also 
considered a highly nonlinear initial value reconnection problem, the Orszag-Tang vortex, 
at these Lundquist numbers (Orszag and Tang 1979, Dahlburg and Picone 1989). Re- 
connection occurs within less than one Alfven transit time, with the attendent formation 
of electric current sheets and vortex quadrupoles, in a manner similar to previously re- 
ported numerical simulations of this configuration. Hence we also conclude that nonlinear 
magnetic reconnection can occur at these Lundquist numbers. We have also run parts 
of the magnetic arcade simulation reported in the present paper with a smaller viscosity, 
M = 30, and have found no significant qualitative change in the results. Let us emphasize 
that these Lundquist numbers, although large by the standards of numerical simulations, 
are still many orders of magnitude smaller than the expected solar values; hence, the 
simulations axe bound to be much more diffusive than the real corona. 

Equations 1 are solved in a three-dimensional Cartesian geometry. We employ the 
following numerical algorithm to discretize the equations: A Fourier pseudospectral dis- 
cretization is used in x a nd y, the periodic directions. Chebyshev collocation is imple- 
mented in the z direction, with the mechanical pressure evaluated on a grid staggered in z. 
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The temporal discretization n a low-storage, third-order, Runge-Kutta - Crank-Nicolson 
scheme, with a backward-Euler pressure correction at each time-level (Zang and Hussaini 
1986). For a grid of 32 2 Fourier modes and 32 Chebyshev modes the code requires about 
2 seconds per time step on the NRL CRAY X-MP, and approximately .8 million words of 
memory. More details on the algorithm aro given in the appendix. 

b) The Initial and Boundary Conditions 

Since the code requires periodicity in the x and y directions we pick as the initial field 
a periodic array of current-free arcades. The initial vector potential is given by: 

A x = A x — 0, (2o) 

A y = cos kx exp (—k(z + 1)) , (2b) 

and the computational domain consists of 0 < x < 2tt,0 < y < 27T, — 1 < z < 1. This i3 
the same potential field studied in 2.5 dimensions by Mikic et al., 1988 who considered 
the evolution of this field under the influence of photospheric shear flows. The initial field 
has translational symmetry along the y direction; however, this symmetry is broken im- 
mediately by the flow imposed at the bottom boundary of the computational box, which 
is supposed to correspond to the photosphere. Unlike previous calculations (Mikic et al., 
1989, Van Ballegooijen 1988) which focused on equilibrium states, we concentrate on coro- 
nal dynamics and do not attempt to obtain an equilibrium or steady-state. Our simulation 
consists of two distinct periods, a “driving” part during which free energy is built up in 
the coronal magnetic field by footpoint motions, and a “decay” part during which the field 
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relaxes back to a current-free state. The field is first twisted with a prescribed footpoint 
motion for 10 characteristic time scales r, where r is the average Alfven crossing time for 
the whole system. Then the footpoint positions axe fixed, and the system is allowed to 
evolve for another 60r. During the driving period, t < 10 r, the horizontal velocities at the 
bottom wall z = —1 axe assumed to be, 


where ip is given by, 



and 


Vy 


tty 

dx ’ 


i p — C sin 2 2 kx sin 3 ky, 


( 3 ) 

( 4 ) 


and C is a constant. Of course during the decay period, t > lOr, v s = v« = 0. Note 
that the form (3) for the footpoint flow corresponds to an incompressible twisting motion. 
The bottom wall is also assumed to be impenetrable, so that v x — 0 throughout the 
simulation. The boundary conditions at the other walls are straightforward. At the sides 
of the computational box we must assume periodic boundary conditions. The top wall, 
z— 1, is assumed to be a rigid, impenetrable surface, consequently v = 0 there. 


The only other quantity that requires boundary conditions is the vector potential A. 
Again periodic conditions are imposed at the sides. The top and bottom conditions on A 
are selected so that the evolution at the boundary corresponds to ideal flow. However, this 
determines only the horizontal components of A, 

^ = v,B" ( 5 «) 
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dAy 

at 


VxJ^X) 


(56) 


where we have used the fact that v z vanishes on the boundary. In principle, no boundary 
conditions need be imposed on A z , but since the numerical code requires them, we simply 
solve the partial differential equation for A z at the boundary (c/. Mikic et a/., 1988). 

!T = (vxB),- ij,. (5c) 

We have also tried extrapolating its ■< r alue from the interior, but the results are not sensitive 
to this assumption. 

Due to the combination of the periodicity restriction and the particular problem under 
consideration, there is an unavoidable redundancy in the calculations. There are two 
arcades in the computational box, one full arcade with its neutral line at the center of 
the box, x = 27r/6, and one half arcade on either side. Throughout the simulation we 
found that both arcades behave identically, (as they should); hence, we will only show the 
central arcade in the following figures. This redundancy will not be present in calculations 
with more realistic magnetic fields and photospheric flows. For example, the photospheric 
flow pattern is perhaps best represented as a random, time evolving flow-field. Such a 
photospheric flow pattern would break the symmetry implied by equations 2 and 4. 

The value of C in equation (4) for the twisting motion is chosen so that the maximum 
velocity at the photosphere is approximately 1/10 of the coronal Alfven speed. Of course, 
this ratio for the driving speed to the characteristic velocity is much greater than in the 
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solar case. Typical photospheric velocities are of order 1 km/s and coronal Alfven velocities 
~ 1,000 km/s so that the actual ratio should be of order 10 -3 . However, it is not possible 
to use such a low ratio in the numerical simulations because diffusion is so much stronger 
in the simulation model than on the reeil Sun. If the driving speed is too low, then the 
effects of the footpoint motions will diffuse away before any twist can be built up in the 
coronal magnetic field; the field is able to “slip” back faster than it is driven. Also, since 
the time step in the code is set by the characteristic velocity, it is simply too expensive 
computationally to use a very low ratio. A value of ~ 1/10 is about the smelliest that can 
be used with our present code. This implies that our results during the driving period are 
likely to be too dynamic compared to the solar corona; however, the simulation should be 
accurate for the decay phase. 

In Figure la the str?am lines of the photospheric driving flow, equivalent to con- 
stant contours of Vs are shown, and in Figure lb, contours of the velocity magnitude, 
|v| = [(dip/dx) 2 + (drp/dy) 2 ] \ The flow pattern has been chosen so that the velocity 
vanishes at the boundary between adjacent arcades and at the neutral line of each ar- 
cade. This ensures that each arcade retains its identity; there is no mixing of flux between 
different arcades or between opposite polarity regions. From Figure 1 we note that the 
velocities are quite concentrated and are oppositely directed on either side of the neutral 
line so that they correspond to twisting a fairly localized flux region. Again due to the pe- 
riodicity requirement there is a redundancy in the velocity pattern, consequently we have 
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two regions of twisted flux. They behaved identically throughout the simulation, which 
lends confidence in the numerical code. 

The footpoint displacements that result from the photospheric velocities are shown in 
Figure 2, where we plot the contours of constant B z on the bottom wall at the beginning 
and the end of the driving period, t — lOr. Note that from equation (2) the initial 
photospheric flux B z = —ks'mkx so that the initial contours are straight lines parallel to 
the y axis (Figure 2a). From the ideal induction equation one finds that B z is convected 
with the flow for an incompressible horizontal motion, 

+ Vpfc • V-B* = 0, (6) 

where v p h is the horizontal velocity at the photosphere. Hence, by following contours of 
B z one can track the footpoint positions. 

There axe several important features of Figure 2. It clearly demonstrates the point 
that for ideal motions the topology of the coronal field can be determined from the bound- 
ary conditions at the photosphere. The fact that the coronal field has been twisted is 
obvious from the photospheric contours. Even if the footpoint motions were much more 
complicated, involving winding or braiding patterns and compressible flow, the topology 
would be recoverable from the contours of the Euler potentials for B at the photosphere 
(Antiochos 1987,1990). It is also evident from Figure 2 that due to the footpoint motions, 
the B z contours have been “stretched out”, and become more concentrated. We conclude 
that the gradients in the footpoint displacements increase with time and, therefore, the 
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scales of the magne .ic structures in the corona must decrease with continued twisting (e.g. 
Van Ballegooijen 1988, Mikic et al., 1989). Note that this occurs even if the scales of the 
photospheric velocities remain constant. This result implies that a photospheric velocity 
with large scales such as that given by (4), will eventually produce structure on the photo- 
sphere, (and consequently in the corona), with scales below the numerical resolution. The 
actual magnitude of the twisting is evident in Figure 2. Each side of the arcade has been 
twisted through a maximum angle ~ 180°, so that the total twist is approximately one 
full rotation, which is about the maximum twist that we can resolve with a 32 3 grid. In 
order to resolve a more complicated motion such as a braiding pattern which requires at 
least three rotations (Antiochos 1990), we would need a much larger grid. 

III. RESULTS OF THE SIMULATION 

The evolution of the magnetic field lines during the simulation are shown in Figures 
3a - 3d (side view), and 4a - 4d (top view). We plot two sets of field lines at four times 
during the run, t = 0, lOr, 20r, and 70r. These field lines were chosen by the position of 
one of their footpoints on the photosphere. We selected two rows of 17 points lying on the 
photosphere, and then traced the field line emanating from these points. The two rows 
are parallel to the y axis and pass approximately through the center of the twisted region 
of the arcade. In the Figures, these points are to the left of the photospheric neutral line 
which lies at the center of the grid, (#,■ = 2 : 17 ). Field lines with footpoints on the inner 
row are shown as dotted lines, those on the outer axe solid. 
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Figures 3a and 4a show the initial current-free magnetic field (equation 2). The 
translational symmetry along the y direction is apparent. The field clearly has a very 
simple topology with all of the field lines from the inner row of footpoints lying well under 
the lines from the outer row. Figures 3b and 4b shows the field at the end of the driving 
phase, t = lOr, when it has the maximum twist. It should be noted that these are not the 
same field lines as in 3a and 4a. The field lines of Figures 3a and 4a have moved from their 
original positions so that their leftside footpoints no longer lie on our two straight rows (see 
figure 2). However, there is nothing special about the particular field lines in 3a and 4a, 
and we find it easier to see the effect of the motions by keeping a fixed, two-straight-rows 
pattern for the leftside footpoints in all the figures. 

Comparing 3a, 4a and 3b, 4b it is evident that the field topology becomes much more 
complex as a result of the twisting. The field lines now wrap around each other so that 
there is no clear distinction between an outer and an inner set on the right hand side of the 
photospheric neutral line. The most striking feature of the twisted field lines is that they are 
all greatly expanded upward from the current-free positions. This result is expected from 
studies of force-free fields (e.g., Aly 1984; Yang, Sturrock and Antiochos 1986). The field 
acts to minimize any stress due to footpoint motions by expanding outward in the corona. 
Note that this effect is not incorporated in models with an initially uniform field confined 
between two horizontal plates (e.g., Parker 1972, Mikic et al, 1989). For a given footpoint 
motion, such models will yield fields with finer current structures than in a more realistic 
open coronal geometry. In fact, our results also underestimate the outward expansion. 
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The boundary conditions that we assume at the top and bottom of an impenetrable wall, 
and the restriction to incompressible flows, equation (lb), all act to suppress an outward 
expansion. 

The expansion plays a crucial role in determining the field’s structure and dynamics. 
We note from Figures 3b and 4b that all the field lines are concave down in the corona. 
There are no “dips” suitable for supporting cool prominence material. If the field lines 
had not expanded upward, we believe that dips would have formed. The reason for the 
lack of dips is that outward expansion tends to decrease the pitch of field lines as they 
twist around each other. It seems likely that more twisting will eventually produce dips, 
on the other hand, Aly (1984) has argued that continual stressing of the field should lead 
to an open arcade. If so then even infinite twisting will not produce dips in the corona. 
Unfortunately, our simulation is not capable of accurately treating more than about one 
rotation so that we cannot address this interesting question at present. 

The topology of the field lines at later times is shown in Figures 3c - 3d and 4c - 
4d. Since the footpoints are held fixed during this phase of the evolution, then these 
field lines would correspond exactly to those in 3b, 4b if the evolution were ideal. It is 
evident from the Figures that the evolution is far from ideal. The right-side footpoints are 
clearly changing with time even though we allow no flow on the boundary. The reason for 
this is resistive diffusion, which permits the field lines to “slip” back to their current-free 
positions. It is important to emphasize that this slippage is not a boundary effect; it is 
occurring throughout the corona. On the contrary, at the boundary the code is able to 
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ensure an ideal, perfectly rigid plate; but in the interior volume the code must use a finite 
resistivity. Note that this slipping is not unphysical; it also occurs in the real corona, but 
at a much slower rate. 

We find that the evolution throughout the decay phase is one of simple diffusion. 
It is evident from Figures 3 and 4 that all the field lines tend to untwist. There is no 
evidence for any type of instability, ideal or resistive. This can be seen ^rom plots of the 
global quantities. We show in Figure 5 the evolution during the whole simulation of the 
logarithm of the kinetic energy (E v = l/2(|v| 2 )), the magnetic energy (Eb = 1/2(|B| 2 )), 
and the mean square electric current ( J = l/2^|j| 2 )). The forcing at the wall causes E v to 
build rapidly to a saturation level. The kinetic energy decreases rapidly, due to the high 
viscosity, at the end of the driving phase (f = 10)and remains small thereafter. Growth of 
E v , typically seen when jets form as a result of reconnection (e.g. Dahlburg, Dahlburg,and 
Mariska 1988), is not observed. The magnetic energy, Eb, also grows quickly during the 
forcing phase to about 1.5 times its initial value. After the forcing is cut off at t = 10, Eb 
exhibits a smooth monotonic decay to a level below its initial value. This reflects the fact 
that the magnetic field relaxes to a different equilibrium, consistent with the new footpoint 
positions exhibited in Figure 2b. In the same way, J increases during the forcing phase, 
as the forcing creates smaller spatial scale structures in the magnetic field. During the 
damping phase, J decays to near zero, indicating the relaxation of the magnetic field to a 
nearly potential state. 
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The nature of the magnetic field evolution can be described in terms of an effective 


magnetic Lundquist number, given by 


S(t) = 


2<|j| 2 ) 

£<W). 


m 


This is also shown in Figure 5. We obtain S(t ) by considering the magnetic induction 
equation in the limit of negligible velocity, taking the scalar product of this equation 
with B, and integrating the resulting magnetic energy density equation over the system. 
Solving for the Lundquist number gives the desired result. Figure 5d shows that S(t) 
returns rapidly to the predicted Ohmic rate, viz., S(t) = 500, when the forcing stops at 
T — lOr. The “jiggling” in the data is due to small errors in computing the denominator 
of S(t), which is postprocessed. Since the magnetic energy fails to decay at a higher rate, 
magnetic reconnection must be absent. 


Another important result of our simulation is that it shows no evidence for current 
sheet formation. Hence, our results are consistant with the studies of force-free equilibria 
(Sakurai 1979; Van Ballegooijen 1988; Mikic et ai, 1989). The lack of current sheet 
formation can be seen from the evolution of J shown in Figure 5. Not only does the 
integrated current exhibit a monotonic decay, but we have also examined in detail the 
spectrum of the current density, and have found that the spectrum decays as would be 
expected for pure resistive diffusion. 
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rv. DISCUSSION 


The results presented in this paper are of a negative character, i.e.,, we do not see 
magnetic reconnection, kinking, or the formation of concave-up magnetic field lines suitable 
for prominence formation. Hence, it appears that, contrary to many models, photospheric 
twisting of a single arcade does not lead to the type of processes required to explain solar 
activity. The key question is whether this conclusion is valid in general, or whether it is due 
solely to the limitations of this particular simulation. The numerical simulation itself has 
three major shortcomings. First, it is too diffusive, especially when compared to the real 
corona. There is little that can be done about this problem. Larger Lundquist numbers can 
be obtained by going to larger grids, but there is no hope at present for obtaining numbers 
significantly larger than ~ 1000. Second, the numerical resolution is barely adequate. At 
the time of maximum twist, lOr, the current spectrum shows a decrease of only 4 orders of 
magnitude over the wavenumber range. Ideally, we would like to have at least a 6 orders of 
magnitude decrease. Since we have used only a 32 3 grid, and much larger grids are possible 
on present computers, we see no difficulty in obtaining substantially better resolution in 
future simulations. A third problem is the incompressibility assumption. This is not valid 
for the solar corona; but, we believe it has not had a significant effect on the results. If 
anything, a compressible simulation should exhibit an even larger upward expansion and 
less of a tendency to instability. 
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There are also limitations with the assumed form for the photospheric motions. We 
believe that the magnitude of the twisting is sufficient. It is evident from Figure 5b that 
the magnetic free energy is 50% the energy of the potential field. Also it can be seen 
from Figures 3b and 4b that the field lines are strongly twisted. Coronal loops rarely 
appear to be so distorted. Although the magnitude of the motions is acceptable, their 
size scale is limited. Due to the difficulty with obtaining adequate numerical resolution, 
we have considered only a large-scale twisting in which the width of the twisting motion 
is approximately equal to the width of the arcade. This situation would correspond to a 
large-scale shearing of a coronal arcade, and is the easiest to resolve numerically. However, 
the photosphere also contains small-scale motions such as granular motions, for example. It 
may be that the small-scale motions play a crucial role in processes such as the formation 
of prominence dips. Simulations with larger grids are needed in order to address this 
possibility. 

Even though our simulation is clearly limited, we believe that at least some of the 
results will hold for general footpoint motions and Lundquist numbers. In particular, 
we expect that for a dipolar magnetic arcade, resistive instabilities will not occur due 
to the effect of line-tying, irrespective of the footpoint motions. This argues against the 
reconnection scenarios proposed by a number of authors. We also expect that, contrary 
to the non-equilibrium model, current sheets will not form. Of course the scale of the 
currents in the corona is determined by the scale of the footpoint displacements, so that 
if the displacements have a very fine scale, then fine scale currents will result. However, 
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our simulation and all the equilibria calculations (Sakurai 1979; Van Ballegooijen 1988; 
Mildc et al., 1989) find that for well-resolved footpoint displacements smooth current 
distributions result. 

The question of ideal instabilities and the formation of concave-up field lines is less 
clear. We find that both are suppressed by the outward expansion of the field. Such 
an expansion will occur for any form of footpoint motions, but it is possible that strong, 
highly-localised twisting can lead to the formation of concave-up field lines and to kink-like 
instabilities in small regions of the arcade. If so, then small scale motions may have im- 
portant consequences for prominence formation and coronal dynamics. This issue requires 
further study. 

The general conclusion of our results is that photospheric stressing of a magnetic ar- 
cade is incapable of producing phenomena such as prominence eruption or coronal heating. 
Some critical ingredient is missing from this model. There are several possibilities. One is 
that emerging flux plays the dominant role in activating the coronal magnetic field. Sim- 
ulation which allow for flux emergence and submergence need to be considered. Another 
possibility is that a complex magnetic topology is required. We have argued (Antiochos 
1989,1990) that a pure dipolar field (i.e., one due to a single positive and a single negative 
polarity region on the photosphere) is somewhat pathological, since it is the only field with 
a well-behaved connectivity everywhere. In a more realistic field which is due to a number 
of different polarity regions at the photosphere, the connectivity must be discontinuous 
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and null points must occur in the corona. In this case current sheets and reconnection will 
readily occur as a result of footpoint motions. We are currently pursuing these possibilities. 

In summary, we believe that the work reported here constitutes only a first step in 
examining the 3D dynamics of coronal magnetic fields. Much more interesting simulations 
await to be done. 


APPENDIX 

COMPUTATIONAL PROCEDURE 

We modify the algorithm of Zang and Hussaini (1986) to include magnetic field terms. 
The governing equations 1-5 arc solved within the computational range ft by means of 
the following time-step splitting: 

v* =v*xw*+j*xB* + if y2y ( in 

Aj=v*xB*-^VxVxA (in ft) 

o 

v* = g* (on 3ft) 

A* = h* (on 3ft) 

v*(x,t„) = v(x,t„) (in ft) 
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A *(x,i„) = A (x,f„) 

(in fi) 

v?* = - VP** 

l 

(in f2) 

V • v** = 0 

(in Ji2) 

v** • n = g • n 

(on 

v** • f = -VP** • f 

(on 5f2) 

v**(x,f*)=v*(x,i*) 

(in 0) 


The first time-level, represented by single-starred variables, is a convective, dissipative 
step. The second level, represented by double-starred variables, is a pressure step which 
enforces the solenoidality of the flow field. 

A staggered grid is employed using the following collocation points: x, = 2ni/N x for 
i — 0, 1,JV* — 1, yj = 2irj/N y for j = 0, l,N y — 1, z* = cos(7r k/N g ) for k = 0,1, N Z) 
and, zif+i /2 = cos[7r(fc + 1/2 )/N z ] for k = 0, 1 , N g — 1 . Velocities and vector potentials are 
defined at the nodes (x{,yj,Zk). Pressures are defined at the shifted points (z,j, yj,z fc+1 / 2 ). 
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PHOTOSPHERIC FLOW FIELD 








NORMAL MAGNETIC FIELD AT PHOTOSPHERE 



Fig. 2 — One quadrant of the normal magnetic field in the photosphere. The same quadrant is shown in Fig. 1. (a) 
t = Or, initial conditions, (b) t = 10r, conclusion of the twisting phase. The figure shows how much the magnetic field is 
twisted from its initial configuration, approximately one full rotation. After t = 10r, the photospheric B z retains this form. 







(c) (d) 


Fig. 3 — One quadrant of magnetic field lines, side view. The same quadrant is shown as in Fig. 1. (a) t = Or, initial 
conditions, (b) t = 10r, conclusion of the twisting phase, (c) t * 20r, partial damping, (d) t * 70r, conclusion of the 
simulations. 
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(C) 


(d) 


Fig. 4 — One quadrant of magnetic field lines, top view. The same quadrant is shown as in Fig. 1. (a) t « Or. initial 
conditions, (b) t - lOr, conclusion of the twisting phase, (c) t - 20r. partial damping, (d) / * 70r. conclusion of the 
simulations. 
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